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I. INTRODUCTION 



Understanding the rheological properties of particle dispersions has been an important 
problem in many fields of science and engineering. When a dispersion is subjected to shear 
flow, the flow properties of the dispersion show a variety of non-Newtonian behaviors such 
as shear thinning and shear thickening. These non- Newtonian behaviors are associated with 
the changing microstructures of the dispersion, and several different physical mechanisms 
for these peculiar behaviors have been proposed. 

In recent years, several numerical methods have been developed to accurately simulate 
particle dispersions, and they are all based on a similar approach, which involves resolving 
the fluid motion simultaneously with the particle motion. We refer to this approach as 
direct numerical simulation (DNS). Recently, we have developed a numerical method, which 
we call the smoothed profile method (SPM), for the DNS of particulate flows.-- In the 
SPM, the Navier-Stokes equation for the fluid motion is discretized on a fixed grid, and 
the Newton's and Euler's equations for the particle motion are solved simultaneously with 
the fluid motion. One simple technique to impose shear flow with the DNS approach that 
maintains conventional cubic periodic boundary conditions is to apply a spatially periodic 
external force to generate a periodic flow profile. We have successfully used a zigzag flow 
profile to impose both steady and oscillatory shear flows in the DNS of spherical particle 
dispersions.- 1 ^ 

When a zero-wavevector shear flow is required, the usual cubic periodic boundary con- 
ditions must be modified to be compatible with a time-dependent shear deformation of the 
simulation cell. Such a modification was proposed by Lees and Edwards^ and is commonly 
used in various simulation studies. The Lees-Edwards boundary conditions can be very 
easily implemented for particle-based simulations such as molecular dynamics simulations. 
However, care must be taken to implement these conditions in continuum grid-based simula- 
tions such as computational fluid dynamics or time-dependent Ginzburg-Landau equations. 
The most useful implementation of the Lees-Edwards periodic boundary conditions for grid- 
based simulations is to solve the dynamic equations in deformed (oblique) coordinates.-~- 
Onuki proposed a general methodology to examine the phase transition dynamics and rhe- 
ology in the presence of shear flow,- and it has been successfully used in several simulation 
studies and particularly for polymeric fluids in shear flow.- — 
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The aim of this short paper is to propose a method to implement the Lees-Edwards 
periodic boundary conditions to simulate dispersions of solid particles in host fluids by the 
combinatory use of the SPM and the oblique coordinates. 

II. METHOD 

In the SPM, the boundary between the solid particles and the solvent is replaced with a 
continuous interface by assuming a smoothed profile. This simple modification enables us 
to calculate hydrodynamic interactions both efficiently and accurately without neglecting 
many-body interactions. The equation governing the dynamics of particle dispersion is a 
modified Navier-Stokes equation: 

[du 1 

P \~dt + ^ U ' V ^ I = V ' a + P ^ p ~ Kp ( Ux ~ ^ 
with the condition of incompressibility V • u = 0, where p is the solvent density, 

a- = -pi + f) f {V«+ (V«) r } (2) 

is the Newtonian stress tensor with a solvent viscosity of rji, and u(r,t) and p(r,t) are 
the velocity and pressure of the dispersion, respectively. A smoothed profile function < 
4>( r ,t) < 1 distinguishes between the fluid and particle domains as well as yields 0=1 
in the particle domain and = in the fluid domain. These domains are separated by 
thin interstitial regions with thicknesses characterized by £. The dispersion density p is 
represented as 

P = C 1 - <f>)Pt + <1>Pp ( 3 ) 

where pf and p p are the solvent and particle densities, respectively. Only neutral buoyancy 
dispersions with p = pf = p p are simulated in the present study. The body force <pf p is 
introduced so that the total velocity field u of the dispersion satisfies u(r) = (1 — </>)uf(r) + 
<f)u p (r), where uj is the fluid velocity and u p represents the rigid motions of the particles. 
The incompressible condition V • u thus ensures V0 • (u p — uf) because both Uf and 
Up satisfy incompressible conditions. The gradient of <p is proportional to the surface- 
normal vector and has a support on the interfacial domains. Therefore, the body force 
<f)f p introduced to satisfy the rigidity of the particles ensure the appropriate impermeability 
boundary conditions at the fluid-particle interface, while the non-slip boundary conditions 
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are imposed automatically by the viscous stress term in the Navier-Stokes equation. More 
detailed explanations and the mathematical expressions for <f> and </>/ p were also detailed in 
our previous papers.-^ 

The last term in Eq. ([T]) represents the external force needed to maintain linear shear flow: 

u x = jy (4) 

where 7 is the shear rate, and K is a constant that determines the amplitude of the external 
force. Here we impose only the zero-wavevector shear flow so that the averaged fluid velocity 
becomes compatible with Eq.(j4]). 

The motion of the z-th particle in a dispersion is governed by Newton's and Euler's 
equations of motion: 

= + *T + a". |-. = « (o) 

= +9? (6) 

where 7*3, v^, and a;, are the position, translational velocity, and rotational velocity of the 
colloidal particles, respectively. M, and J, are the mass and the moment of inertia, and 
//* and nf are the hydrodynamic force and torque exerted by the solvent on the colloidal 
particles, respectively^, gj and gf are the random force and torque, respectively, due 
to thermal fluctuations. The temperature of the system is defined such that the long- 
term diffusive motion of the colloidal particles reproduces the Stokes-Einstein rule.— >^ ff 
represents the potential force due to direct inter-particle interactions such as through the 
Coulombic and Lennard- Jones potentials. 

Eqs.fll]), (G3), and (EJ) are solved simultaneously in the SPM. However, this task is not easy 
with an ordinary periodic boundary condition because Eq.([T]) depends explicitly on y, which 
leads to a violation of the translational invariance. This problem can be eliminated by using 
oblique coordinates. Fig. [1] represents a schematic illustration of the present coordinate 
transform. At a time t = t , a spherical solid particle is located in a solvent in Fig. [1] 
(a) where the solvent is discretized into square grids in an ordinary rectangular coordinate 
system. In Fig. [1] (b), the grids are deformed due to the shear flow that is applied for 
t > to while the shape of the solid particle is unchanged. The same situation is depicted in 
a transformed (oblique) frame in Fig. [1] (c) where the grid has not moved (i.e., it remains 
rectangular), but the shape of the solid particle changes over time due to the shear flow. 

4 



To formulate the oblique coordinate transformation based on tensor analysis, we began 
by redefining the covariant basis Ei and contravariant basis E i in oblique coordinates rather 
than using the expressions shown in the literature.- - - Fig. [2] provides a definition of the basis 
vectors. Using a rectangular unit vector, Ei and E l are expressed as 

Ei = e x E 1 = e x - ^te y 

E 2 = A fte x + e y E 2 = e y (7) 
E 3 = e z E 3 = e z 

where e a is the unit vector in the a(= x, y, z) direction in the original rectangular coordinate 
system. We can obtain contravariant (covariant) vector components A 1 (Aj) using A ■ E % 
(A ■ Ei). The positional vector r = xe x + ye y + ze z is transformed from the rectangular 
coordinate expression r to the oblique coordinate expression r as follows: 

r = xe x + ye y + ze z 

= (r ■ E v )Ei + (r ■ E 2 )E 2 + (r • E 3 )E 3 (8) 

= x l Ei + x 2 E 2 + x 3 E 3 = r, 

where the contravariant components (x l ,x 2 ,x 3 ) are expressed as 



x 1 = x — 'yty 

* 2 

x — y 



(9) 



x = z 



and the time in oblique coordinates is expressed as i — t. Each contravariant component 

is transformed to a covariant component by using the metric tensors Gij = Ei ■ Ej and 
G l i = E l ■ E^ . Then, the transformation can be expressed as 

A i = G ij Aj (10) 

Ai = Gi 3 & (11) 

The physical quantities in Eq. are transformed as indicated below: 

p(r,t)=p(r,t) (12) 

0(r,t) = 0(r,t) (13) 

u(r, t) = u(r, t) - 7?/e x (14) 
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4>f P (r,t) = <t>f p (r,t) (15) 

In the oblique coordinate system, u satisfies the standard periodic boundary conditions while 
u satisfies the Lees-Edwards periodic boundary conditions in the rectangular coordinate 
system. The contravariant components (u l ,u 2 ,u 3 ) of u are expressed as 

u 1 = u x - jtUy - 72/ 

u 2 = u„ (16) 



u 3 = u, 



where (u x , u y , u z ) are the rectangular components of u. The contravariant components of 
f p are (0/p, 4>f 2 , 4>f 3 ) and can be expressed as 

- j i 



- - jt<t>fy 

fv = <t>f v P 
fp 3 = 0/ P 2 



(17) 



where (f£ , fg, /*) are the rectangular components of / p . 

The differential operators in oblique coordinates are defined by 



^ dx 1 ^ dx 2 + ^ dx 3 



d d . d 

5=«+™& (19) 



Therefore, the Laplacian operator in oblique coordinates is expressed as 



dx 1 J \dx 2 dx 1 J \dx 3 r 
Using these formula, Eqs. (JTJ and (j2J) are rewritten in oblique coordinates as 

Hlif 

I ^ J (21) 

= V ■ <r + p4>f p -pju 2 Ei - Kp(u l + 7« 2 )^i 



and 

**(r, f) = -G*p(r , t) + w {g-^ + G-J^} (22) 

with the incompressibility condition V • u = 0. Because Eq. ([21]) and u satisfy the periodic 
boundary conditions in all directions, a fast Fourier transformation (FFT) can be safely used 
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to solve the Poisson equation, which is needed to determine p with the incompressibility 
condition. In Appendix 1, detailed explanations are given on how to solve Eq. (]2ip with 
the incompressibility condition in the oblique coordinate system with using the spectral 
(Fourier) method. 

When 7 = jt — 1, the positions r = (x l , x 2 , x 3 ) on an oblique grid with 7 can be mapped 
onto the identical positions r = (x, y, z) on the original rectangular grid with 7 = using the 
operation x = x 1 + ct 2 , y = x 2 , z = x 3 . The shear strain 7 is then reset to 0.- Repeating this 
process allows us to perform stable numerical calculations over a long period with keeping 
< 7 < 1. The above coordinate transformation based on the tensor analysis leads to the 
same expression for the Laplacian A as that of a previous study- However, a difference 
arises between the differential operators V for which our formal transformation derives a 
much simpler expression as shown in Eq. (fTKj) . 

We calculate the dynamics of solid particle dispersions in shear flow by following these 
steps: 

i) The fluid velocity field in the oblique coordinate system at a new time t = nh is 
calculated by integrating Eq. (j2"Tl) over time with 0/ p = as 



u* = u n ~ 1 + 

t„_!+h r . /1 \ i (23) 

ds 



V ■ fid- - iiuj - {K{u l + 7M 2 ) + 27M 2 } Ei 

while satisfying the incompressibility condition V • u* = 0. Here, the superscript n denotes 
the time step, and h is the time increment. In Appendix 1, detailed explanations are given 
also on how to solve Eq. (!2B1 with using the spectral method. 

ii) The velocity field u* is transformed into rectangular coordinates u* using the inverse 
transformation expressed as 
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u x = u 1 + 7M 2 + 7X 
u y = u 2 (24) 
u z = u 3 . 

iii) The motions of colloidal particles are only calculated in rectangular coordinates. The 
position of each colloidal particle is calculated by 

rtn-l+h 

r? = r^~ 1 + / v^ds (25) 



iv) Using the momentum conservation between colloidal particles and the solvent, the hy- 
drodynamic force and torque acting on each colloidal particle are computed with volume 
integrals within the particle domain as 

cH _ P I ( „.* „.n-l 

and 



f? = % I dr[fi («• - u;- 1 )] (26) 



nf = P - J dr[(r - r,) x # (u* - v^ 1 )] (27) 

where , Up~ 1 (r) = ^j</>"(f) (vf + x (r — 7*j)J is the correct velocity field within 
the particle domain in which <p ~ 1. The space integrals in Eqs. fl26l) and fl27j) are carried out 
by summations over grid points in actual computations, however, there occur grid mismatch 
between u* which is supported on oblique grid points and other variables (</>" and 
^Wp" 1 ) which are supported on rectangular grid points r^j^. We determine values of u* on 
rectangular grid points r^j^ by linear interpolation as described in detail in Appendix 2. The 
translational velocity and rotational velocity of each colloidal particle are then calculated as 

i rtn-i+h 

jr\ (t? + fr + 9?)ds (28) 



and 

= + /r 1 / («f + tf) <** (29) 

v) To ensure the rigidity of the particles and the appropriate non-slip boundary conditions 
at the fluid/particle interface, the body force <pf v is calculated as 

(u™ - u*) 1 
0/p = ^ ^ - -Vpp. (30) 

The correcting pressure p p is determined to make the resultant total velocity incompressible. 
This leads to the Poisson equation of p p : 

Ap P = P ^ (31) 

We then transform 0/ p into oblique coordinates (f>f p using Eqs. (TT5l) and ( 1T7T) . 

vi) Finally, we obtain the correct fluid velocity field as: 

ii n = ii* + 0/ p /i. (32) 



Repetition of steps i) through vi) provides a complete procedure to perform the DNS of 
colloidal dispersions under shear flow. 

We can calculate the stress tensor of the dispersion (s) and the dispersion viscosity 
V = ( s xy)/i in the following manner where (•••) denotes averaging over space and time. 
The equation governing the dispersion is formally written as: 

^(pu) = V • <r dis - Kp(u x - 7?/K (33) 

By comparing Eq. ([T|) with Eq. (133|) . we get the formula 

V-<r dis = V-<t + P 0/ p . (34) 

The full stress tensor s of the flowing dispersion is then defined by introducing a convective 
momentum-flux tensor explicitly as 

s = <r dis - puu. (35) 

The definitions of cr dls and s are identical to the definitions in our previous paper.- Now, 
we can evaluate the average stress tensor of the dispersion (s) directly from Eqs. ( 134"|) . (|35|) . 
and 5cr = s — er as 

<«> = 



with the volume V = L x L y L z where Li is the system size in i-direction. (•••)* denotes time 
averaging over steady state. In the derivation of the second formula, we used a second rank 
identity. If we substitute Eq. (134|) into the third formula, then we obtain the fourth formula. 
The fifth formula can be obtained by assuming that the system is in a steady state in which 

(I (P»)>« = (I + u ■ V (^)>* = and (I (P«)>* = o. 
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(36) 



III. RESULTS 

Using the method described above, we calculated the high- and low-shear limiting vis- 
cosities of colloidal dispersions for various volume fractions of particles $. The particles 
interact via a truncated Mie (m,n) potential with m = 36 and n = 18.— 



where r is the distance between the centers of a pair of particles. The parameter e character- 
izes the strength of the interactions, and a represents the diameter of the colloidal particles. 
The lattice spacing Sx is taken to be the unit of length. The unit of time is given by pf5x 2 /r] 
where rj = 1 and pf = p p = 1. The system size is L x x L y x L z = 64 x 64 x 64. Other 
parameters are set as follows: a — 8, £ = 2, e = 1, r\ — 1, Mj = na 3 /6, and h = 0.067. The 
temperature is ksT = 7. The range of shear rate is 1.0 x 10~ 4 < 7 < 0.1. 

The inset of Fig. [3] shows the dependence of the Newtonian viscosity on the volume 
fraction $ when $ < 1. The present simulation data show very good agreement with Ein- 
stein's viscosity law. Fig. [3] shows the dependence of the low-shear limiting viscosity (closed 
symbols) and the high-shear limiting viscosity (open symbols) on the volume fraction. Our 
simulation data for both high- and low-shear limiting viscosities show good agreement with 
the experimental results of van der Werff et al— Previously, Brady theoretically predicted 
the behavior of the low-shear limiting viscosity— Our simulation data show good agreement 
with Brady's prediction over a wide range of volume fractions < $ < 0.55. Ladd analyzed 
the behavior of the high-shear limiting viscosity using Stokesian dynamics.— Our simula- 
tion data agree well with Ladd's simulation data and also with the theoretical results of 
Beenakker.— 

Finally, we add some comments on the differences between the present method using 
Lees-Edwards boundary condition and the previously proposed method using zigzag velocity 
profile.— We simulated a single spherical particle in shear flow using the two methods without 
thermal fluctuation. The volume fraction is 0.001. Figure shows the ratio of angular 
velocity a; of a spherical particle to the applied shear rate 7 as a function of 7. Although the 
data using the present method tend to be slightly smaller than the data using the zigzag flow, 
deviations of both data from the analytical value w/7 = 0.5 remain small within numerical 
errors of the methods. Figure [6] shows the intrinsic viscosity [r/] of the dilute dispersion as 




(37) 



10 



a function of shear rate 7. The simulation data using the present method almost perfectly 
follow onto the Einstein's prediction [rj] = 2.5, while the data using zigzag flow slightly 
overestimate [77] because of unphysical kinks of the zigzag flow profile. This problem is not 
very serious when the shape of dispersed particles is spherical and the size of the particle 
is much smaller than the distance between two kinks. Serious problems, however, occur if 
this method is applied to non-dilute dispersions of chains or rods, for example. The present 
method using Lees-Edwards boundary condition is free from this problem. 

IV. CONCLUSION 

We presented a generic methodology for performing DNS of particle dispersions in a 
shear flow using oblique coordinates and periodic boundary conditions. The validity of 
the method was confirmed by comparing the present numerical results with experimental 
viscosity data for particle dispersions over a wide range of the parameters $ and 7 that 
include nonlinear shear-thinning regimes. An important advantage of the DNS approach 
over other approaches such as Stokesian dynamics is its applicability to particle dispersions 
in complex fluids. In fact, electrophoresis of charged colloids^ and particle dispersions in 
nematic liquid crystals^ have already been calculated using SPM. Our methodology can 
also be applied to simulate particle dispersions in viscoelastic fluids simply by replacing the 
Newtonian constitutive equation to more complex ones such as Maxwell model. 
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APPENDIX 1 

In this section, we describe how to solve Eq. (12 lj) with the incompressibility condition 
in an oblique coordinate system using Fourier spectral methods. The Fourier and inverse 
Fourier transforms are defined as 
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A(r) 



A(k) = J A{r) exp(— ifc ■ r)dr 
1 



A(k) exp(ik ■ r)dk, 



(38) 
(39) 



where k is the wavevector of the oblique coordinate system. In k space, we can express the 
spatial covariant derivative as 



dA(r) 
dx a 



— > ik a A(k). 



(40) 



where k a is a covariant component of k. 

Using these relations, we modify Eq. ( 123]) from r space to k space. This equation is 
solved in k space. First, by substituting Eq. fl22|) into Eq. f l23|) . we obtain the explicit 
equation represented by 



u (r) —u(r)+ 

^-i+ h r A „ ) 

(«(r) ■ V)w(r) - V^-^ + uAu(r) - {Kiu 1 ^) + -fu 2 (r)) + 2^u 2 (r)} E x 



t n -l 



P 



(41) 
ds. 



where v = rjf/p is the kinetic viscosity. Using a Fourier transform, the form of Eq. (14T1) in 
A; space is written as 



u*(k) = u(k)+ 

"tn-l+h 



F(k) - uk 2 u(k) - \K{u l (k) + 7« 2 (fc)) + 2^u 2 {k)\ E x 



ds, 



(42) 



where 



F(k) = / (ii(r) ■ V)-u,(f) exp(— ifc • r)dr. 



The bracket 



(43) 

A(k) {= A(k) ■ (^1 — ||JJ denotes taking the orthogonal part to A; and 
this operation corresponds to imposing the incompressibility condition V • u(r) = (or 
equivalently k ■ u(k) = 0). Because p(r) is automatically determined by imposing the 
condition of incompressibility, we can safely neglect this term. 

Using the method describe above, we calculate u*(f) from u(r) as shown in Eq. ff23l . 
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APPENDIX 2 



An arbitrary position vector in the oblique coordinate system is defined as 

r lu = (iEi + jE 2 + kE 3 )5x (44) 

with arbitrary integer numbers i,j, k, while a position vector in the rectangular coordinate 
system is defined as 

ri,j,k = (ie x + je x + ke x )Sx (45) 

with integer numbers i,j,k, where Sx represents the lattice spacing. In general, the two 
position vectors ~. ^ and are not compatible with each other. To transform from ^ 
to r ijifc using Eq. i must equal i — jj. However, jj is not always an integer since 7 is 
defined between and 1. We thus perform interpolation of the variables to overcome this 
problem. 

Fig. |4] shows a lattice discordance between the rectangular and oblique coordinate sys- 
tems, is the location vector in rectangular coordinates, and and rjj-g are the 
location vectors in oblique coordinates. Using liner interpolation, we estimate the veloc- 
ity field in rectangular coordinates u(r it j^) from the velocity field in oblique coordinates 
u ifi_ijk) and it(fjj^). From the liner interpolation, u(r*jj )fc ) along the straight line is 
given by the equation 

= ^fri 77\ u (n-i,u) + • ., —\ u ( r i,u)- ( 46 ) 

I i,j,k i—l,j,k\ I i,j,k i— l,j,k\ 

When using liner interpolation, artificial diffusion may arise. To check the reliability 
of the present method, we calculate the angular velocity u and intrinsic viscosity [77] = 
lim$_ 5 . ( ? 7 — Vf)/® f° r a dilute dispersion of spherical particle for which analytical solutions 
are available. As already shown in Figs. [5] and |6l the present simulation data agree very 
well with analytical solutions indicating that the effects of artificial numerical diffusion are 
not serious. 
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Rectangular Oblique 

FIG. 1. A schematic illustration of the present coordinate transformation. In (a), a spherical solid 
particle is in a solvent, which is discretized into grids in an ordinary rectangular coordinate system, 
at a time t = to- Since the shear flow is applied for t > to, the solvent (grids) is convected by the 
flow while the shape of the solid particle is unchanged. Such a situation is depicted in the original 
(experimental) frame in (b) and also in a transformed (oblique) frame in (c). The transformation 
between (b) and (c) is defined by Eq. ([7|) . 
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0.1 0.2 0.3 0.4 0.5 

volume fraction <j> 

FIG. 3. The behavior of the viscosity rj as a function of the volume fraction <J>. The open symbols 
represent the high-shear limiting viscosity, and the closed symbols represent the low-shear limiting 
viscosity. The open and closed circles correspond to our simulation data, whereas the triangles 
correspond to experimental results^ The solid line is Brady's theoretical prediction^ and the 
dotted line is a fitting curve obtained from previous experimental and simulation^ data. The 
inset indicates a comparison of our present simulation data with Einstein's viscosity law (dashed- 
line) in the small volume fraction regime <E> < 0.01 where the viscosity exhibits simple Newtonian 
behavior. 
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FIG. 4. A schematic illustration of the lattice discordance between the oblique and rectangular 
coordinates. is the location vector in a rectangular coordinate system. rj —1 •■ ^ and f j ~- £ are 

the location vectors in the oblique coordinate system. 
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FIG. 5. The behavior of the ratio of angular velocity oj to the shear rate j as a function of 7. 
Open circles indicate the results of the previous method. Closed circles indicate the results of the 
present method. The solid line corresponds to the analytical solution. 
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FIG. 6. The behavior of the intrinsic viscosity as a function of shear rate 7. Open circles indicate 
the results of the previous method. Closed circles indicate the results of the present method. The 
solid line corresponds to Einstein's viscosity law. 
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